function H_tp1_peers = mean_peers(H1,H2,PS1,PS2, gamma,fp,set_kids2,fp_parfor,t) 

%if fp_parfor.j~=3
%scale_util = fp.scale_util(1) + fp.scale_util(2)*fp_parfor.n_simulation/100 + fp.scale_util(3)*( (fp_parfor.n_simulation /100)^(2));
%scale_util = fp.scale_util(fp_parfor.j);
%else
%scale_util = 0;    
%end

U_links1 = gamma(1).*ones(size(H1))  +  ...
           gamma(2).*log(H1) + gamma(3).*log(H2) + gamma(4).*(log(H1) - log(H2)).^(2) + ...
           gamma(5).*( (log(H1) - log(H2)).^(2) ).*(log(H1) - log(H2)>0).*PS1;                                                 
       
U_links2 = gamma(1).*ones(size(H1))  +  ...
           gamma(2).*log(H2) + gamma(3).*log(H1) + gamma(4).*(log(H1) - log(H2)).^(2) + ...
           gamma(5).*( (log(H1) - log(H2)).^(2) ).*(log(H2) - log(H1)>=0).*PS2;         
               
Prob_Links_temp = (1./(1+1./exp(U_links1))).*(1./(1+1./exp(U_links2))) ;
    
%Links = binornd(1,repmat(Prob_Links,[1,1,fp.tot_draws_mcintegration_network]),[size(Prob_Links,1),size(Prob_Links,2),fp.tot_draws_mcintegration_network]) ;      

Prob_Links = repmat( Prob_Links_temp, [1 , 1 , fp.tot_draws_mcintegration_network]) ; 

Links = ( Prob_Links >= fp_parfor.peers_backward(:,:,:,t));        

Child_util = repmat( U_links1, [1 , 1 , fp.tot_draws_mcintegration_network]) ; 

Child_util = Child_util.*Links;

peers = repmat(set_kids2',[size(Prob_Links,1) ,1, fp.tot_draws_mcintegration_network]) ;
         
temp = Links.*peers;

H_prime = nansum(temp,2)./sum(temp~=0 &  isnan(temp)==0,2) ;  

H_prime = reshape(H_prime,[size(H_prime,1),size(H_prime,3)] ) ;                        

H_prime_sq = H_prime.^(2);       
        
H_tp1_peers.H_prime = H_prime;
H_tp1_peers.H_prime_sq = H_prime_sq;
H_tp1_peers.U_links1   = Child_util ;